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The population dynamics of single- and bi-exciton states in semiconductor nanocrystals is modeled 
numerically in the presence of Coulomb coupling between single- and two-exciton states and a 
dissipation channel in order to study the transient bi-exciton population that occurs in an optically 
excited semiconductor nanocrystal. The results show that the system evolution strongly changes 
if the dissipation is included. In a certain range of parameters, the growth of the exciton number 
(MEG process) is fast (on picosecond time scale) and the following decay (Auger process) is much 
slower (hundreds of picoseconds). In some cases, the maximum occupation of the bi-exciton state 
increases when dissipation is included. The dynamics of an ensemble of nanocrystals with a certain 
size dispersion is studied by averaging over the energy of the bi-exciton state which can be different 
for each single nanocrystal. The validity of Markov and secular approximation is also verified. 



I. INTRODUCTION 

One of the possible ways of improving the efficiency of 
the existing solar cells is to exploit the process of multiple 
exciton generation (MEG) in semiconductor nanocrystals 
(NCs)W Such an effect consists in generation of two or 
more electron-hole pairs by a single high energy photon 
and thus converts the excess above-bandgap energy into 
useful current. This process is enabled by Coulomb cou- 
pling between single-pair (exciton, X) states and two-pair 
(biexciton, BX) states in a NC (or, in general, between 
states with n and n + 1 pairs) and consists in an intra- 
band relaxation of a carrier (typically an electron, due to 
larger energy scales of confined states in the conduction 
band) accompanied by a creation of a new electron-hole 
pair (an inverse Auger process). In this way, the excess 
energy obtained by an electron upon absorbing a high- 
energy photon is not dissipated in a phonon relaxation 
processes and becomes available for photovoltaic conver- 
sion. 

The initial experimental results, showing very high 
values of the quantum efficiency of photoconversion in 
various systemsP^, were subsequently reinterpreted^*^ 
based on the growing understanding of the experimental 
difficulties that might lead to overestimating the achieved 
numbers of excitons per single absorbed photorf^^. 
Nonetheless, m ore re cent experiments on real NC-based 
solar cell devicesSsEH do provide a direct proof of the use- 
fulness of this process in solar energy conversion. The- 
oretically, the description of the X and BX spectrum 
and the X-BX couplings that are essential for the MEG 
process h as be en proposed using the methods of density 
functional 20 -* 2 ^ pseudopotentiaP 2 ^ 2 ^, tight binding^"22l, 
and k-p theory®^. 

Along with the investigation of these structural prop- 
erties, much attention has been devoted to the carrier 
dynamics in a nanocrystal under optical excitation at 
energies high enough to generate multiple excitons, in 
particular to the role of decoherence and relaxation. 
These studie s included dynamical simulations of few-level 
modelpSlMI as we \\ as f many-level models aiming at re- 



producing the density of X a nd BX states in the high- 
energy sector of a nanocrysta l 29 * 31 * 3 -^^. In many cases, 
dissipative effects are included in these models on a phe- 
nomenological level and expressed by a number (usually 
small) of dephasing rate d 26 * 30 * 31 * 33 !. In this way, the mul- 
tiple exciton generation could be describ ed as a process 
competing with exciton relaxatio n 1 30 * 31 * and, in certain 
cases, suppressed by coupling to the dissipative environ- 
ment (typically considered to be phonons^. 

In this paper, we study the time evolution of the X- 
BX system within a minimal, three-level model that ac- 
counts both for the impact ionization and Auger recom- 
bination in the presence of dissipation. Expressing the 
couplings to the environment in a generic form in terms 
of a physically motivated set of spectral densities allows 
us to characterize the emerging coupling to the Coulomb- 
correlated X-BX eigenstates and to discuss the depen- 
dence of the rates of various phonon-assisted processes 
(relaxation and impact ionization) on the Coulomb cou- 
pling itself. We show that, on the general level, the dissi- 
pative MEG process is determined by the same couplings 
to the dissipative environment as the carrier relaxation 
and dephasing. Furthermore, we find out that the system 
dynamics realizes various dynamical scenarios, depend- 
ing on the alignment of the X and BX levels and on the 
relation between the level spacing and the spectral prop- 
erties of the coupling to the environment (in particular, 
the high-frequency cut-off of the spectral density). As 
we show, the presence of dissipation considerably mod- 
ifies the system dynamics and, in many cases, increases 
the efficiency of the MEG process. We study also the 
role of the excitation conditions (pulsed, continuous wave 
or incoherent thermal) and show that the strong differ- 
ences between the system kinetics under different excita- 
tion conditions^ are washed out by dissipation. Finally, 
we assess the validity of Markov and secular approxima- 
tions for the description of dissipation-assisted MEG in 
nanocrystals. 

The paper is organized as follows. In Sec. [TTJ we de- 
fine the model (Sec. II A), des cribe the Master equations 
for the system evolution (Sec. II B), and discuss the for- 
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FIG. 1. (Color online) Relative energy of BX states in vicinity 
of selected bright X states with quantum numbers as indicated 
as a function of the NC radius around Ro = 3 nm. The 
numbers on the right denote degeneracy. 



mal structure of the carrier-environment coupling in a 
NC (Sec. |II C ). In Sec. Ill the results of our simulations 
are discussed: first the dynamics in the Markov limit is 
studied (Sec. Ill A I and then non-Markovian corrections 
are discussed (Sec. IIIB). Finally, Sec. IV concludes the 
paper. 



II. MODEL 

A. The system 

Although the density of BX states in a nanocrystal 
is very high, only a few tens of them are coupled to a 
given bright X stated Moreover, usually only a few 
coupled BX states lie in the vicinity of the X state while 
vast majority of them is relatively distant. In Fig. [T] 
we show two selected examples of the spectral positions 
of BX states relative to selected X states to which they 
are coupled (the quantum numbers for the selected ex- 
citon states are the same for the electron and the hole 
and are shown in the figure). This approximate result 
is obtained within the single-band envelope function ap- 
proximation with Coulomb interactions included in the 
lowest order for an In As nanocrystal with the radius R 
close to Ro = 3 nrr i 37 l 38 [ As can be seen, only one or a 
few BX states appear in the spectral vicinity of a given X 
state. Thus, the essential features of the MEG dynamics 
can be expected to be determined by impact ionization 
within groups of a few states. 

Consider an optically active (bright) excited X state 
of a NC. At lower energies, other X states are present 
to which the carriers can relax. For our purpose, it is 
sufficient to consider one such level. We will consider 
situations in which a BX state Coulomb-coupled to the 
bright X state is present in between the two X states, 
which is a common situation for highly excited X states, 
where the coupled BX states are rather dense. Taking 
into account the selection rules for intcrband Coulomb 
couplings^, it is rather unlikely that this BX state will 
also be coupled to the other, lower X state. Therefore, for 
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FIG. 2. (Color online) (a) The schematic diagram of sys- 
tem states and couplings. The green arrow shows the lase 
excitation while the black arrows represent the environment- 
induced transitions. The Coulomb coupling between X and 
BX configurations is shown with a red arrow, (b) The diagram 
of system eigenstates with the optical excitation paths (green) 
and environment-induced transitions (blue), (c) Graphical 
representation of the processes involved. 



our dynamical modeling, we consider a four-level model 
of a nanocrystal, as shown in the Fig. [2|a). Here the 
state G refers to the ground state (empty nanocrystal), 
A and B denote the two X states and 2 is the BX state. 
We assume that the BX state 2 is coupled only to the 
state B by a Coulomb coupling V which is taken to be a 
real parameter. 

The Hamiltonian of the carrier system is then 



H = e A \A)(A\ 



£b\B)(B\ 
VV{\B){2\ 



-e 2 |2)(2| 



(1) 



where €a,£b and e 2 are the energies of the states A, B 
and 2 respectively. We set the energy of the ground state 
to be zero. Moreover, the system is excited by a clas- 
sical light pulse which has the frequency close to the 
G < — > B transition. By standard selection rules, this 
pulse couples the ground state only to the X states. We 
assume that only the state B is bright. The correspond- 
ing Hamiltonian of the excitation is 



#ias = \f{t) {\G)(B\e 



iQt 



\B)(G\e 



-int\ 



(2) 



In addition, the system undergoes dissipative dynamics 
due to the interaction with its environment. We do not 
assume any specific form of the coupling to the environ- 
ment and aim at a model which is independent of the 
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nature of this coupling. The only restriction on this cou- 
pling is that in the absence of Coulomb-induced mixing 
between the X and BX configurations no interband pro- 
cesses are allowed [see Fig. [2][a)]. Thus, the interaction 
with the environment can be described by the following 
Hamiltonian 

H int = \A)(A\K AA + \B)(B\K BB + |2)(2|^ 22 

+ \A)(B\K AB + \B)(A\K B a, (3) 

where lZ a /3 (with a, (3 — A, B, 2) are certain environment 
operators with the property lZ a p = TVa a - 

The system evolution is described in the basis of eigen- 
states of H : \G), \A), |+) and |-), where |+) and |-) 
result from the Coulumb coupling between the X state 
\B) and the BX state |2), 

|+> = cos(0/2)|B) +sin(0/2)|2), 
|-)=-sin(0/2)|B)+coB(0/2)|2), (4) 

Here, 9 is the mixing angle, defined by tan 9 = 2V/(cb — 
€2), which is close to for weakly mixed X and BX states 
and goes to ±7r/2 if a nearly degenerate pair of X and 
BX states is strongly coupled. The energies of these 
states are E± = E ± U, where E = (es + £2) /2 and 
U = [(cb — £2) /4 + V 2 ] 1 / 2 . In the eigenstate basis, the 
interaction Hamiltonian can be written as 

Hint = (TijTlij, (5) 

i,3=A,± 

where er.y = and 

9 9 
K A + = cos - Kab, K A - = - sin - TZ A b, (6a) 

£±± = \ (Kbb + TIaa) ± 2 cos 9 (Kbb - TIaa) , (6b) 

n+- = ^sw9{K 2 2-Kbb), (6c) 

with TZ± = Kji. 

B. Evolution: Master equation 

While Markov approximation has commonly been used 
to model the dissipation effect on the MEG process in 
nanocrystals (oft en on the level of phenomenological de- 
phasing ratesJPsHSHMEIID its validity is not obvious for 
the present system. Indeed, as observed in experiment 
and reproduced by our simulations discussed below, only 
the Auger recombination phase of the system dynamics is 
slow, while the initial impact ionization dynamics takes 
place on much shorter, picosecond time scales. Therefore, 
in this paper, we will compare the Markovian dynamics 
with a more general approach, where the reservoir mem- 
ory is taken into account. In both cases, the evolution 
of the system will be described in the density matrix 



formalism by solving the appropriate quantum Master 
equation. 

Thus, the starting point for our modeling of the system 
evolution is the time-convolutionless Master equation for 
the reduced density matrix of the charge subsystem in 
the NC in the lowest order. In the interaction picture, 
this equation has the forrrPH 

1 /■* 

P = Tr R dr [H iDt (t), [H int (r) , p(t) ® p R }} , (7) 

where Hi nt (t) and p are the interaction Hamiltonian 
and the reduced density matrix of the nanocrystal in the 
interaction picture, pr is the density matrix of the reser- 
voir at thermal equilibrium, and Tr R denotes the partial 
trace over the reservoir degrees of freedom. Upon sub- 
stituting the interaction Hamiltonian from Eq. |5]) and 
taking the partial trace this yields 

~(f\ _ _ i(wi-Wj+Wk-Wi)t 



*Fijkl{t) WijVkiP~{t) ~ cr fei/ 5(i)cr.y] + h.c, (8) 
with Ui = Ei/h and the time-dependent rates 




where IZij(s) denotes the operator TZij in the interaction 
picture with respect to the reservoir Hamiltonian and 
we have neglected the imaginary parts of the rates that 
describe reservoir-induced energy shifts. 

The reservoir correlation function ("memory func- 
tion") (Hij(s)lZki) is related to the spectral density 

Rijufa) = ^ J dte^KytyKu), (10) 

for k,l = A, ±, which fully characterizes the dissipa- 
tive coupling to the environment. In the same way, spec- 
tral densities TZ a p^s{uj) are defined in terms of correlation 
functions between the operators R a p in the original basis 
[Eq. Q3p] . If the reservoir correlations decay on a certain 
time scale (reservoir memory time) then, on longer time 
scales, the rates become constant and equal to 

Fijki(t) Tijki(oo) = 2iTRij k i {ui - uj k ) . 

Moreover, as follows from Eq. Q, in the absence of de- 
generacy, the rates other than Tijji and Tujj contain an 
oscillating factor and can be assumed to be small if the 
separation of the levels is large and the overall system dy- 
namics is slow. Therefore, it is common to use the secular 
approximation where the terms containing such oscillat- 
ing rates are neglected. The rates of the type ^ajj at 
long times are proportional to the corresponding spec- 
tral density at zero frequency, which vanishes in many 
common situations (super-Ohmic reservoirs and Ohmic 
reservoirs at zero temperature, see below). Thus, one 
reaches the commonly used Markov approximation with 
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the evolution equation in the Lindblad forirPH, which in 
the Schrodinger picture can be written as 

p = --[H , p] + ^2T ij (a 3l pa l3 - -{(Tycr^, p}),(ll) 



where I\j = r^^oo) 



C. Environment: Spectral densities 

Obviously, the spectral densities defining the transi- 
tion rates are related to those in the original basis. In 
particular, 



1 



R 



A±±A 



( W ) 



sin 2 6»[i?2222(w) 

— R22Bb{u) — -RbB22(w) + RbBBb(u)\, 
± COS 6) R abb A, (12) 



with Ujiij(u>) = TZijji(uj). 

It is clear that the relevant spectral densities Rijji de- 
scribing the transitions between Coulomb-correlated X- 
BX configurations are combinations of the spectral den- 
sities Rapjs that describe dephasing and interaband re- 
laxation between X and BX states. Interestingly, the 
transition between the two Coulomb-mixed states, that 
is, the impact ionization process, described by Rz^±±zf: is 
entirely related to the diagonal couplings between the 
original states and the environment. Obviously, the 
corresponding rate vanishes for small X-BX mixing as 
2 ~ [V/(es — £2)] 2 - On the other hand, the transition 
to the other X state A (relaxation or impact ionization) 
is governed by the off-diagonal couplings, which are re- 
lated to intraband relaxation between these states. If the 
diagonal and off-diagonal couplings are of similar magni- 
tude (as it is the case, e.g., for carrier-phonon couplings), 
then the impact ionization and Auger recombination are 
formally suppressed by a factor sin 2 6 or sin 2 9/2 as com- 
pared to relaxation. However, as we will show below, 
what really matters is the frequency dependence of the 
reservoir spectral density which can make the impact ion- 
ization process favorable, depending on the X-BX level 
alignment. 

The details of the system-environment coupling may 
vary for different nanocrystal systems and, to our knowl- 
edge no general microscopic theory has been proposed for 
its exact treatment. Hence, in most of our simulation, we 
take the spectral densities in the original basis in the sim- 
plest Ohmic form, R a /3 7 s(u>) — a a p 7 s [^b( w ) + 11 J( u)), 

for 



IIIB 



where J{uj) = ( Jquj/wq) exp[— (lu/luq) 2 ] (see Sec 
comparison with a super-ohmic case) 

Here, jib(w) is the Boltzmann distribution function, 
Jo is the overall magnitude of the dissipation and ujq 
is the cut-off frequency. In the simplest case of lowest- 
order acoustic phonon processes, ujq is on the order of 
R/c, where c is the speed of sound, which yields a 



value on the order of ps _1 . This can be different if 
multiple-phonon processes are included. The values of 
the coefficients a a p 7 g follow from the physical nature 
of the carrier-environment coupling: The spectral den- 
sities i?Qct/3^(w) result from diagonal couplings between 
the carriers in the NC and their environment. It seems 
reasonable to assume that these couplings are at least 
approximately proportional to the charge density (this 
is true, e.g., for carrier-phonon couplings as well as for 
Coulomb couplings to fluctuating background or impu- 
rity charge). Thus, typically, IZ22 = 27Z-BB and, in con- 
sequence, a 2 2BB = a B B22 = 2, a 22 22 = 4 and = 1 
otherwise. 

As we will show, the system dynamics in the dissipa- 
tive MEG process depends to some extent on the exci- 
tation conditions. Nevertheless, instead of including the 
electromagnetic field explicitly in our simulation we note 
that the excitation may fall essentially in three classes: a 
short (spectrally broad) pulse excites an optically active 
(bright) X state, a long (spectrally narrow) pulse excites 
selectively an eigenstate of the system and broad band 
thermal radiation excites an incoherent mixture of sys- 
tem eigenstates, proportionally to their overlap with the 
bright X state (see Appendix for details). Thus, as the 
initial state we take the state \B), corresponding to an 
ultra fast coherent excitation (a broad band laser pulse) , 
a pure |+) or |— ) state for a narrow band laser pulse or 
a mixture of the eigenstates |+) and |— ), corresponding 
to incoherent excitation by thermal radiation. 



III. RESULTS 
A. Dissipative MEG dynamics 

Several parameters play a crucial role in the dissipa- 
tive MEG process: the energy differences between the 
X state A, the excited state B and the BX state 2, the 
Coulomb coupling between the X and BX states, and the 
parameters governing the dissipative relaxation process 
(the magnitude Jq and the frequency cut-off u>q). In this 
section, we study the dynamics of the dissipative MEG 
process for various energetical alignment of the states, 
assuming fixed values of the dissipation parameters, and 
compared the results to the case without dissipation. 

In the absence of Coulomb coupling, there is no mech- 
anism for a transition to the BX state (this state is com- 
pletely decoupled), hence, only relaxation between the 
states B and A takes place. This is shown in Fig. [3] (the 
dynamics in this case is independent of the excitation 
conditions). In this case, the BX state is indeed never 
occupied and the total number of excitons, N x , remains 
equal to one. The only occurring process is the intraband 
occupation transfer from the initial state B to the dark 
X state A. 

The variation of the state occupations and the average 
number of excitons in the presence of a realistic Coulomb 
coupling (Vb = 1 meV) is shown in Fig. [4] for the case 
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FIG. 3. (Color online) The variation of the state occupations 
Na, Nb and N2 and the average number of excitons N x in 
the absence of Coulomb coupling at T = K. 



of coherent ultrafast excitation. As mentioned above, 
under these excitation conditions, the initial system state 
is \B), hence iV x = 1. When the bright state \B) and 
the BX state |2) are close to each other, Fig.Hja,b), the 
Coulomb interaction leads to strong mixing of the states 
\B) and |2). As a result, transition between the resulting 
eigenstates are efficient which means that the dissipative 
impact ionization is very fast as manifested by the rapid 
growth (below 1 ps) of the occupation of the state |2) 
and the corresponding increase of the average number of 
excitons (the lower eigenstate is still predominantly bi- 
excitonic). Since the initial state \B) is a superposition of 
system eigenstates, there is an intense oscillation at the 
beginning of the process but its amplitude goes to zero 
in about 10 ps. Later on the system relaxes to the state 
\A), which corresponds to the Auger recombination. The 
BX state decays in this process on a time scale of about 
150 ps. The slow rate of the Auger process is due to 
the relatively large energy distance to the state A, which 
makes the relaxation to this state ineffective. 

For another configuration, Fig. Efc,d), when the BX 
state is shifted closer to the state \A), the relaxation be- 
havior does not change significantly. However, the max- 
imum occupation of the state |2) is much lower. This 
results from the larger energy spacing between the states 
which is beyond the frequency cut-off of the spectral den- 
sity, which considerably suppresses the relaxation to the 
state |— ) which is a predominantly BX state. Thus, in 
the competition between the impact ionization (transi- 
tion to the state |— ) or, almost equivalently, to the BX 
state 1 2)) and the usual relaxation (transition to \A)), the 
latter starts to dominate. In addition, the similar energy 
spacings £5 — 62 and £2 — £A makes the Auger recombina- 
tion from the state 2 to the state A much more effective 
in comparison to the impact ionization. In this case, the 
oscillation amplitude decays in a slightly longer time of 
a few tens of picoseconds. In both of these configura- 
tions, the impact ionization competes strongly also with 
the Auger recombination, the latter being faster in the 
second configuration. 

In the absence of dissipation, J = 0, for both the 
above mentioned configurations, the average number of 
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FIG. 4. (Color online) The occupations of the system states 
(right panels) and the average number of excitons (left panels) 
for Vb — 1 meV for 2 possible alignments of the energy levels 
at T = K. In the left panels, the average number of excitons 
in the case without dissipation. Coherent ultrafast excitation 
is assumed here. 
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FIG. 5. (Color online) As in Fig.Hbut at T = 50 K (a,b) and 
T = 300 K (c,d). 



excitons oscillates constantly (green lines in Fig. |4ja,c)). 
However, as one could expect, the amplitude of the os- 
cillation is different for various configurations due to dif- 
ferent degrees of mixing between the eigenstates. Note 
that for the first alignment (Fig. |4|a,b)), the number of 
excitons generated in the presence of dissipation is larger 
than that achieved without dissipation (taking the aver- 
age of oscillation in both cases). 



6 




20 40 60 80 100 120 140 20 40 60 80 100 120140 

Time (ps) Time (ps) 

FIG. 6. (Color online) The dynamics of (a) the average num- 
ber of excitons (b) the system state occupations for incoherent 
excitation by thermal radiation at T = K. In (a), the aver- 
age number of excitons in the case without dissipation is also 
shown (green line). 

At higher temperature, T — 50 K, Fig. |5ja,b), the 
relaxation dynamics does not change considerably but, 
the oscillation amplitude decays in a shorter time. At 
T = 300 K, Fig. |5j;c,d), the initial average number of 
excitons is slightly decreased and the oscillations vanish 
quickly. The final average number of excitons goes up at 
higher temperatures because of nonzero thermal occupa- 
tions. 

Under incoherent excitation, according to the Fermi 
golden rule, the system eigenstates are excited propor- 
tionally to their coupling to the light, that is, to the ad- 
mixture of the bright state \B). Since in the presence of 
the Coulomb coupling, the eigenstates |+) and |— ) are su- 
perpositions of the original X state \B) and BX state |2), 
they are characterized by the average number of exciton 
between 1 and 2. Hence, the initial number of excitons 
may exceed 1. The system dynamics in this case, shown 
in Fig. [6] for one of the level alignments, is very similar to 
that following a coherent excitation but no oscillations 
are seen since no coherence between the eigenstates is 
present. There is again a strong competition between 
the impact ionization, the X relaxation and the Auger 
recombination. As a result, the occupation of the BX 
state again grows rapidly and then decays completely on 
a longer time scale. 

Under these excitation conditions, in the absence of 
dissipation, the average number of excitons would remain 
constant but in the presence of dissipation it increases 
considerably to a maximum value and then decays. It 
is clear that the maximum average number of excitons 
resulting from the dissipative MEG dynamics in this case 
exceeds that following the excitation in the absence of 
dissipation. This is due to the dissipative transition to 
the predominantly biexcitonic state |— ), which develop 
on a time scale much shorter than the subsequent Auger 
transition to the lowest state \A). 

It can be seen by comparing Fig. [6] and Fig. [4] that for 
incoherent excitation, the value of the average number 
of excitons is equal to the mean value of the oscillations 
that could be seen in Fig. [4] The same holds true also 
for the other level alignments, not shown in Fig. [6j 
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FIG. 7. (Color online) The variation of state occupation and 
the average number of excitons as a function of time at T = 
K. (a,b) for the initial state |+); (c,d) for the initial state 
|— ) (narrow band excitation). In (b,d), the initial lOps of the 
evolution is expanded. 



Under the third possible excitation conditions, when 
a spectrally narrow light field is used, a light beam is 
tuned to excite just the state |— ) or |+) (see Appendix 
for details). These two cases are compared in Fig. [7j One 
can see that the only difference between starting from |+) 
or |— ) is the behavior of the BX occupation at the very 
beginning of the evolution, which rises during the first 
few ps for the initial state |+) (Fig. |7|b)), which can 
easily be explained. Since the energy of the state |+) is 
higher than the energy of the state |— ), starting with the 
state |+) will cause the transition from |+) to |— ). After 
this initial redistribution of occupation, the occupations 
of the states do not depend on the initial conditions. This 
is in contrast with the case without dissipation where the 
average number of excitons remains about 1.4 and 1.6 in 
the case of starting from |+) and |— ), respectively, in 
accordance with the composition of these eigenstates in 
terms of the X and BX states. Obviously, no oscillations 
can be observed in the evolution of the system which was 
initially prepared in one of the eigenstates |+) or |— ). 

In order to understand the dependence of the ob- 
served dynamics on the energy spacing between the lev- 
els, we will now discuss the case when the energy differ- 
ences between the states are 1.4 times larger than the 
set of parameters in Fig. [4] The dissipation parameters 
J = 1 ps -1 and ujq — 2 ps^ 1 are kept unchanged. The 
results of the simulations are shown in Fig. [8] The im- 
pact ionization now appears in a shorter time interval 
(about 10 ps). If the state \2) lies close to the state \B) 
[Fig. [8ja,b)], because of the large energy difference be- 
tween the states |+) , |— ) and \A) (beyond the cut-off 
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FIG. 8. (Color online) The state occupations (right panels) 
and the average number of excitons (left panels) for 2 align- 
ments of the energy levels in the case of increased energy 
difference between the states (1.4 times larger than Fig.[4| at 
T — K. In the left panels, the upper envelope of the oscil- 
lations of the average number of excitons in the case without 
dissipation is also shown (green dotted line). Coherent ultra- 
fast excitation is assumed here. 



energy huio), the rate of the Auger recombination is de- 
creased as compared to the previous case and the BX 
occupation persists much longer. The oscillation ampli- 
tude is much lower than in the previous case (due to 
a smaller mixing between the X and BX states) and de- 
cays rapidly (in a few ps). As it is depicted in Fig.[8jc,d), 
when the BX state gets closer to the state A, the degree 
of the impact ionization decreases dramatically. Besides, 
oscillation amplitude decays in a much longer time of a 
few hundreds of picoseconds. The achieved number of ex- 
citons is much lower than in all the previously discussed 
cases. 

At higher temperature, Fig. [9j the impact ionization 
process takes place on the time scale of several picosec- 
onds (decreasing with increasing temperature) and is 
then followed by a slow Auger recombination on the order 
of nanoseconds for this set of parameters) . 

In a realistic sample, one deals with an ensemble of 
nanocrystals with a certain size dispersion. To obtain the 
dynamics for the whole nanocrystal ensemble, we average 
our results over the energy of the BX state |2) which 
can be different for each single nanocrystal. We use a 
Gaussian distribution for the value of £2, 



/(e 2 ) = 



1 



2-KCT 



(13) 



1.0- 10" 2.0- 10- 

Timc (ps) 




l.o-icr 2.0-itr 

Time (ps) 



FIG. 9. (Color online) As in Fig. [8] but for T = 50 K (a,b) 
and T = 300 K (c,d). 



keeping tA and es constant. 

The results are shown in Fig. [10] for different sets of 
parameters. In Figs. 10 'a,c,e), we fix the mean BX en- 
ergy e 2 o at 4 meV and the energy es at 4.5 meV above 
tA and show the results for three different values of a. In 



where AE = £2 — ^20 and a is a standard deviation, while 



Figs. 10 b,d,f ) , a small standard deviation, a — I meV is 
chosen and the results for two different level alignments 
are shown. 

Comparison of Fig. 10 with Figs. [4] and [5] shows that a 
small inhomogeneity of the level alignment in the ensem- 
ble (red lines in Figs. 10 'a,c,e)) does not change the over- 
all system dynamics. The amplitude of the oscillations 
at T = is reduced due to ensemble dephasing but the 
average trend is almost exactly the same. This is true for 
all the level alignments studied here, and at all temper- 
atures, as can be seen in Figs. 10'b,d,f). On the other 
hand, increasing the inhomogeneity reduces the ampli- 
tude of the initial peak of the average number of exci- 
tons, in particular at higher temperatures, and leads to a 
nearly featureless time dependence of this quantity after 
the initial ultrafast growth. For larger inhomogeneities, 
also the overall (long time) value of the average number 
of excitons is reduced. This suppression of the MEG effi- 
ciency in the ensemble in our model is due to the fact that 
in a strongly inhomogeneous ensemble the contribution 
from NCs with very distant levels becomes larger. 



B. Non-Mar kovian corrections 

In this section, we discuss some further technical as- 
pects of the modeling of the MEG dynamics in NCs. We 
assess the corrections due to the reservoir memory (be- 
yond the Markov approximation) and study the differ- 
ences between Ohmic and super-Ohmic reservoir models. 
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FIG. 10. (Color online) The variation of the average number 
of excitons as a function of time at three different tempera- 
tures. Left panels: for two different values of a. Right panels: 
for two different energies of the BX state for a — 1 meV. 
(a,b) T = K, (c,d) T = 50 K and (e,f) T = 300 K. Ultra 
fast excitation conditions is assumed here. 



In Fig. 11 we compare the simulation results for one 
selected set of system parameters obtained from the 
Lindblad equation (red solid lines) and from the non- 
Markovian TCL equation (blue dashed lines). In order 
to understand the role of various approximations made 
on the way to the Markovian description in terms of the 
Lindblad equation, we present also results following from 
two kinds of secular approximations to the TCL equation 
(green dotted and gray dash-dotted lines). In the approx- 
imation labeled "SEC.l", only the rates IYyi(t) appear- 
ing in the Lindblad equation are kept (but no Markovian 
approximation is made, which is reflected on the level of 
the TCL equation by the time dependence of the rates). 
In the second secular approximation, denoted "SEC. 2", 
also the other subset of secular rates, Tajj(t) is retained. 
At long times, these rates become proportional to the 
spectral density at zero frequency, hence they tend to 
zero in the zero temperature limit. Therefore, at low 
temperatures, they are only important in the initial pe- 
riod of the dynamics (for times on the order of 1/wq), 
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FIG. 11. Comparison of the simulation results using different 
Markovian and non-Markovian approximations. 



when non-Markovian effects are dominant. 



As one can see in Fig. 11 'a), at zero temperature the 
corrections to the Markovian dynamics are rather small 
and essentially amount to a reduced amplitude of the 
oscillations observed during the first few tens of picosec- 
onds after the excitation. Both the time-averaged value 
of this oscillating exciton number as well as the long-time 
asymptotics are nearly the same in both cases. This is 
also visible in Fig. [TT^b) , where the initial period of time 
is shown in more detail. A close look at the curves re- 
veals that the difference between the Lindblad and TCL 
results is due to initial damping within a few picoseconds 
from the initial time, while the subsequent evolution is 
characterized by the same exponential damping of the 
oscillations and decay of the populations in both cases. 
This initial damping is due to the fact that the short- 
time values of the non-Markovian damping rates involve 
interaction with the whole reservoir spectrum, while in 
the Markov limit only the resonant modes are involved, 
corresponding to the relatively low values of the spectral 
density in its high-energy tails in the present case. 

In Fig. [TTfb), we have also displayed the results ob- 
tained from the non-Markovian TCL equations in which 
the non-secular (oscillating) terms have been removed. It 
is clear that both classes of secular terms must be kept in 
the non-Markovian description even though only one of 
them appears in the Lindblad equation which yields quite 
accurate results. Otherwise, the initial, non-Markovian 
damping is overestimated by an order of magnitude. Still, 
however, the trend is reproduced correctly. Similar effect, 
although with a larger quantitative difference in the oscil- 
lation amplitudes, is seen in a system with lower energy 
spacing between the levels (not shown here). 

A similar comparison for T = 30 K, presented in 



Fig. 11 'c), shows a larger difference between the Marko- 
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FIG. 12. Comparison of the Markovian and non-Markovian 
approximations for an Ohmic and super-Ohmic reservoir. 



vian and non-Markovian modeling results. Now, the 
Markov approximation overestimates the damping of the 
initial oscillations and yields a higher number of excitons 
(the latter depends on the level alignment; we have ob- 
served an opposite situation for smaller inter-level spac- 
ings, not shown here). Much more interesting is the anal- 
ysis of the role of secular vs. Markov approximations 
in this moderate temperature case, shown in Fig. [TTJd) . 
Here, both secular approximations yield the same long- 
time asymptotics but differ in the way the initial oscil- 
lations are reproduced: if only the terms entering in the 
Lindblad equation are kept the oscillations vanish com- 
pletely, as is also the case for the Lindblad description. 
The other, more complete secular approximation, repro- 
duces the oscillations qualitatively correctly, although 
with a shift of the average trend. Even in this case, how- 
ever, the quantitative characteristics of the MEG process 
are very similar in all the approximations. 

So far, in all our simulations we have assumed a reser- 
voir with Ohmic spectral characteristics (J(w) ~ u> at 
u) — > 0). This is the standard choice for generic model- 
ing of the dynamics of an open quantum system in the 
absence of any detailed knowledge about the particular 
reservoir in question and has also been employed in some 
studies of the dissipative MEG process 3 ^. However, some 
specific reservoirs are known to have other spectral char- 
acteristics. In particular, the three-dimensional acous- 
tic phonon reservoir shows a super-Ohmic behavior with 
J(w) ~ lj 3 at u) —> 0. This leads to vanishing spec- 
tral densities R(uj) at uj — and, as mentioned above, 
to suppression of a class of damping rates in the long 
time limit. It seems interesting to study how this quali- 
tative difference affects the system dynamics both in the 
Markov limit and in the non-Markovian model. To this 
end, in Fig. [l2j we show the simulation results for the 
same system as in Fig. [TT] but with a super-Ohmic spec- 
tral density J(uj) = Jq(oj / ojq) 3 exp(— lo 2 /lu'q) with the am- 
plitude J Q and cut-off frequency ujq chosen in such a way 
that the Markovian transition rates from the states |+) 
and |— ) to the state \A) are the same as in the previ- 
ous Ohmic case. The zero temperature results presented 
in Fig. 



12 a) show a very similar behavior to the Ohmic 
case at the same temperature [Fig. [TT|(b)] but much less 



sensitivity to the approximation used. In particular com- 
parison between the results obtained with various secular 
approximations shows that the dynamics is determined 
by the Lindblad rates Tijji . This can be attributed to the 
reduced magnitude of the spectral density in the vicin- 
ity of lj = 0, which reduces the role of the other secular 
rates not only in the long time limit but also already for 
shorter times. This role of the secular terms of the sec- 
ond type is larger at higher temperatures [Fig. (l2^b)], 
where neglecting them leads to strongly overestimated 
damping. Remarkably, the non-Markovian corrections 
are much smaller in the super-Ohmic case and both the 
average trend of the evolution of the exciton number as 
well as the long-time value are nearly the same here, in 
contrast to the Ohmic model (especially at higher tem- 
peratures). 



IV. CONCLUSIONS 

We have formulated and studied a few-level model of 
dissipative multiple exciton generation and relaxation dy- 
namics following a photon absorption in a semiconductor 
nanocrystal. We have accounted for the interaction with 
the reservoir by introducing a physically motivated set of 
spectral densities which allowed us to relate the impact 
ionization and Auger relaxation rates to the diagonal and 
off-diagonal carrier-reservoir couplings, respectively, and 
to highlight the role of the interband Coulomb coupling 
for the magnitudes of the rates governing these two pro- 
cesses. 

With this model, we have solved the Markovian quan- 
tum Master (Lindblad) equation to investigate the im- 
pact of dissipation on the evolution of the single- and 
bi-exciton occupation. We have shown that the system 
evolution strongly changes if the dissipation is included. 
In many cases, the maximum average number of electron- 
hole pairs (i.e., the efficiency of the MEG process) is 
increased if dissipative transitions are allowed and can 
be close to 2. Thus, dissipation can play a construc- 
tive role in the MEG p rocess and is not neccessarily a 
competing factor^ 30 * 31 !. In a certain range of parame- 
ters, the growth of the exciton number (MEG process) is 
very fast (on picosecond time scale) and the following de- 
cay of the biexciton population (Auger process) is much 
slower (on the time scale of hundreds of picoseconds), 
which means that such a dissipative dynamics following 
an ultrafast excitation cannot be excluded based on the 
observed very fast time scale of the process^. In addi- 
tion, the differences between the system dynamics under 
various excitation conditions, which are present in the dy- 
namics of an isolated nanocrystal, are washed out by fast 
reservoir-induced dephasing. We have verified also that 
similar dissipation-related features in the system kinetics 
are observed in an inhomogeneous ensemble of nanocrys- 
tals. 

We have studied also the sensitivity of the modeling re- 
sults to various formal characteristics of the model. We 
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have shown that the dynamics very strongly depends on 
the position of the high-frequency cut-off of the reser- 
voir spectral density. The simulated dynamics depends 
to some extent on the chosen class of the reservoir mod- 
els but the observed differences between the Ohmic and 
super-Ohmic models are mostly of quantitative charac- 
ter and change neither the qualitative features of the dy- 
namics nor the quantitative expectations for an overall 
MEG yield in a nanocrystal ensemble. We have investi- 
gated also the role of non-Markovian corrections to the 
system dynamics. Although the evolution found from 
non-Markovian equations differs in some cases from that 
obtained in the Markov limit these discrepancies mostly 
have the form of oscillations that are present only during 
the first few tens of picoseconds after excitation and are 
not expected to affect the overall quantitative predictions 
for the MEG yield. 



where po — \G)(G\. The occupations of the |+) and 
|— ) states appear in the 2nd order term, assuming a 
Gaussian envelope for /(<), 



/(*) 



2ttt 



(18) 



one finds the density matrix elements corresponding to 
the coherent excitation in the form 



(±|p|±> = 
(+|p|-> = 



1 1± cos ( 



2ttt 2 

-IH+) = - 



1 



2V27TT 



•sin f/e 2 



(a^+a 2 



(19) 
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VI. APPENDIX: THE DENSITY MATRIX 
ELEMENTS CORRESPONDING TO VARIOUS 
EXCITATION CONDITIONS 

In this appendix, we derive the state occupations after 
optical excitation in various excitation conditions. First, 
we consider the interaction between a strong coherent 
field (as a pulsed excitation treated classically) and our 
4 levels system. A single mode radiation source, such a 
laser, will produce an electromagnetic wave with ampli- 
tude £o(t) and frequency Q, 

£(t) =5 (t)cos0f. (14) 

The interaction Hamiltonian between this electromag- 
netic wave and our system is 



-d ■ £{t) = f(t) cosSli(|G)(fl| + h.c), (15) 



where we defined f(t) — —d ■ £o(t). In the eigenbasis 
(|+),|— )) and in the rotating wave approximation, wc 
have 

hnt = h mm- 



2' 
-|GX- 



|e iA+t cos ■ 



le^-'sin- 



h.c. 



(16) 



where A± = Q, — is the detuning from the transition 
energy. The system state after the optical pulse up to 
the 2nd order in £o{t) is 



Po 



dt[H lnt ,po\ 



1 



— oo 
oc 



dt 



dT [flint [flint, A)]], (17) 



where r is the pulse duration. 

In the case of a narrow band excitation condition, when 
A±t ^ 1, all the exponents in Eq. ( [l9| vanish except for 
the one corresponding to A± = 0. Thus the only non- 
zero element will be (— |p|— ) or (+|p|+) corresponding to 
the laser tuned to A_ or A+, respectively. On the other 
hand, for a broad band excitation condition (short pulse) , 
A±r <C 1, all the exponents in Eq. (19 1 are almost equal 
to 1. After inverting Eq. (2J and substituting to Eq. (19) 
one finds (B\p\B) = 1 and~(2|p|2) = 0. Hence, under this 
conditions, only the state \B) is excited. 

Second, let us consider broad band thermal radiation. 
Then an incoherent mixture of system eigenstates is ex- 
cited, so that there is no coherence between the |+) and 
|— ) states. The interaction Hamiltonian is described by 



flint = d ■ £ = ^2gk\(bk\ + 



(20) 



x (\G)(A\ +cos ( - j |G)(+| -sin I - ) \G)(-\ + h.c), 



where bux and b\, x are photon annihilation and creation 
operators respectively. The occupation of the system 
states resulting from this kind of excitation are propor- 
tional to the corresponding transition states, which can 
be found using the Fermi golden rule. Since the states 
|+) and |— ) are very close compared to the photon en- 
ergy, the difference in the photon density of states and 
coupling magnitude is negligible and one finds up to a 
constant. 



(±M±) 



l±cos6» 



<+H-> = HH+) = o. 



(21) 
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